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1 ■ INTRODUCTION  ! 

) 

The  problem  of  determining  the  charge  distribution  on  a 
perfectly-conducting  flat  disc,  and  hence  the  capacitance  of  such 
a disc,  is  a classical  one.  By  dividing  a square  plate  into  square 
subsections.  Maxwell  obtained  a numerical  solution  of  the  Integral 
equation  satisfied  by  the  charge  density.  A similar  method  (i.e. 
the  method  of  square  subareas)  was  used  by  Reltan  and  Higgins  [2] 
and  Harrington  [3]  to  determine  the  capacitance  of  rectangular  plates. 

I 

A numerical  solution  using  a variational  approach  has  been  considered 

by  Noble  [4] . However,  the  method  of  square  subareas  is  obviously  ; 

not  ideal  for  plates  of  arbitrary  shape.  The  present  report  de- 
scribes a numerical  procedure  for  solving  the  Integral  equation  j 

satisfied  by  the  charge  density  on  an  arbitrarily-shaped  disc, 

> 

and  hence  for  determining  its  capacitance.  j 

The  method  employs  a method  of  moments  approach  in  which  the  | 

charge  distribution  on  the  disc  is  approximated  by  pulses  defined  | 

over  quadrilateral  subdomains.  These  subdomains  are  generated  auto-  ^ 

I 

matically  either  by  the  use  of  the  Zienkiewicz-Philllps  isoparametric  ’ 

coordinate  procedure  [5],  or  through  the  use  of  the  transflnite 
blending  transformations  of  Gordon  and  Hall  [6].  An  advantage  of  the  use 
of  quadrilateral  subdomains  is  that  fewer  expansion  functions  are  re- 
quired, compared  to  a method  of  moments  procedure  based  on  triangular 
subdomains.  The  method  described  is  applied  a variety  of  discs. 
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2.  THE  INTEGRAL  EQUATION  AND  ITS  REDUCTION 


The  charge  density  0(x,y)  on  the  disc  ^ satisfies  the  well- 
knovm  singular  Integral  equation 


1_  [f  o(.X|y) 
JJ  D 


where  V is  the  constant  electric  potential  on  the  disc,  e the 

o 

permittivity  of  free  space  and  D - D(x' ,y';  x',y')  the  distance 
between  the  source  point  (x,y)  and  the  field  point  (x',y'),  is 


given  by 


D - [(x’  - x)^  + (y'  - 


By  dividing  the  disc  into  n quadrilateral  subdomains 
» r ■ l,2,...,n,  and  defining  n pulse  expansion  functions 

fj.(x,y)  ■ 1,  (x,y)  e fj.(x,y)  = 0,  (x,y)  t 

r = 1,2, . . . ,n. 


the  unknown  charge  distribution  can  be  represented  approximately 


o(x,y)  - I fj,(x,y)  , I 

r»l 

where  the  are  constants  to  be  determined.  If  we  now  insert  (2) 
in  (1)  and  apply  a point-matching  procedure  to  determine  the 
potentials  at  the  n points  (x'  , y ' ) 6 , s - l,2,...,n,  we 

8 S 8 

obtain  a system  of  n simultaneous  equations: 


* r 

1 l 
! 


i c A - V , 
‘‘1  sr  r 
r-1 


where 


C . 1 [[  dxdy 

sr  4it£  JJ  D(x',  y';  xo 


Once  the  Interaction  matrix  [C  ] Is  determined,  the  constants 

SIT 

are  readily  obtained  by  solving  (3) , so  that  the  charge  density 
a(x,y)  Is  found  from  (2). 

In  order  to  evaluate  the  Integral  In  (4),  the  quadrilateral 
Is  subdivided  Into  four  triangles  having  the  match-point 
(x^  , y^)  e as  a common  vertex,  the  Integral  being  evaluated 


piecewise  so  that 


c y c 

^sr  4Tie  , ^srk  ’ 


o k=l 


where 


dxdy 

. D(x'  , y'  ; x,3 


Thus  the  calculation  of  the  couplings  depends  ultimately  on 

(1)  the  choice  of  the  subdomains  H and  (11)  the  evaluation  of 

r 

Integrals  of  the  type  appearing  In  (5b) . 

If  Is  the  area  of  triangle  then  using  (2),  an  esti- 
mate of  the  capacitance  of  the  disc  Is 

. n 4 

C - ^ ( I I A A ) . 
r-1  k-1  ^ 


3.  AUTOMATIC  SUBDIVISION  PROCEDURES 


3.1.  Introduction 

In  trying  to  subdivide  a domain  for  a method  of  moments  solu- 
tion of  an  electromagnetic  field  problem,  it  is  Important  that  the 
subdivision  process  should  involve  a minimum  amount  of  data  input 
whilst  incorporating  at  the  same  time  a predetermined  flexibility 
both  in  the  choice  of  the  number  of  subdomains,  as  well  as  in  the  vari- 
ation of  the  density  of  the  match  points.  Although  non-automatic 
(or  manual)  subdivision  procedures  can  be  applied  to  any  domain, 
and  in  particular  to  simple  domains  such  as  the  rectangle,  the 
circle,  etc.,  and  the  parameters  describing  the  domains  so  gener- 
ated can  then  be  input  into  the  computer,  in  practice  the  amount  of 
labor  rapidly  becomes  prohibitive  as  the  number  of  the  subdomains 
increases.  Besides,  the  manual  method  becomes,  in  effect,  a trlal- 
and-error  procedure  for  a domain  of  arbitrary  shape.  In  order  to 
overcome  these  problems,  the  idea  of  the  automatic  generation  of 
nodes  and  subdomains  has  engaged  the  attention  of  many  workers,  notably 
those  interested  in  the  use  of  finite  element  methods. 

Research  on  automatic  mesh  generation  has  proceeded  along  three 
main  fronts.  One  of  these  is  the  technique  introduced  by  Fukuda  and 
Suhara  [7]  in  which  the  aim  is  to  produce,  by  a search  procedure,  a 
triangulation  of  the  domain  of  interest  in  which  the  triangles  are 
'well-formed'  (i.e.  they  do  not  possess  angles  that  are  very  acute, 
a triangle  shape  criteria  being  defined  by  a pre-set  parameter).  In 
this  method  the  interior  nodes  are  initially  generated  randomly  and 
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their  positions  are  later  modified  to  exclude  triangles  which  do  not 
meet  the  shape  criterion.  This  method  was  later  improved  by 
Cavendish  [8],  the  aim  being  the  construction  of  triangular  sub- 
domains  whose  shapes  departed  as  little  as  possible  from  the 
equilateral.  In  a very  recent  improvement,  Shaw  and  Pitchen  [9] 
have  succeeded  in  obtaining  a scheme  which  yields  equilateral  tri- 
angles everywhere  in  the  domain  except  for  some  of  the  triangles 
which  have  two  of  their  vertices  on  the  domain  boundary. 

The  second  approach  is  that  presented  by  Zienkiewlcz  and  his 
co-workers  [5,10-13],  which  uses  an  isoparametric  one-to-one  trans- 
formation to  map  a square  onto  a curvilinear  quadrilateral.  This  method 
has  been  extended  by  them  to  multiply-connected  regions,  curved  surfaces 
and  three-dimensional  domains  by  the  use  of  zones. 

Finally,  there  is  the  Gordon-Hall  technique  [6,14,15]  in  which 
the  unit  square  is  mapped  univalently  (i.e.  invertibly)  onto  a given 
region  (together  with  its  boundary).  This  method,  which  like  that  of 
Zienkiewlcz  and  his  school  is  essentially  interpolatory,  uses  transfinite 
blending  functions  so  that  in  mapping  the  unit  square  onto  the  domain  of 
Interest,  the  interpolatory  function  which  extends  to  the  interior  of 
the  domain,  the  function  defining  the  domain  boundary,  matches  this 
latter  function  at  a non-countable  number  of  points  on  the  boundary.  The 
Gordon-Hall  technique  of  generating  curvilinear  subdomains  is,  in  this 
sense,  an  extension  of  the  Zienkiewlcz  method. 

Apart  from  the  above  methods,  there  are  others  which  use  specially- 
defined  natural  coordinates  and  iterative  schemes  [16],  or  Isoparametric 
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blending  transformations  similar  to  those  Introauced  by  Zlenklewlcz  and 
Phillips,  and  Gordon  and  Hall,  but  specialized  to  the  generation  of 
three-dimensional  meshes  [17],  as  well  as  a search  procedure  which 
alms  at  triangulating  plane  polygons  [18]. 

For  the  computation  of  the  capacitance  of  discs  of  arbitrary 
shape,  it  is  Important  to  choose  an  automatic  mesh  generation  pro- 
cedure which  Incorporates  parameters  so  as  to  satisfy  the  following 
criteria: 

(1)  there  should  be  a systematic  method  for  identifying 
and  counting  the  subdomain  and  their  vertices; 

(11)  it  should  be  possible  to  choose  from  the  beginning  the 

number  of  points  at  which  the  charge  density  is  computed, 
since  this  choice  determines  the  order  of  the  interaction 
matrix  and  can  also  be  used  to  estimate  the  capacitance 
resulting  from  the  use  of  an  infinite  number  of  expansion 
functions  (i.e.  the  'asymptotic  capacitance'); 

(ill)  allowance  should  be  made  as  much  as  possible  for  the 

variation  of  the  charge  density  on  the  disc,  such  that 
more  nodes  are  available  near  edges  (where  the  charge 
density  becomes  infinite)  than  elsewhere. 

Both  the  Gordon-Hall  and  Zienklewlcz-Phillips  procedures  provide 
the  simplest  means  (based  on  a Cartesian  system)  for  meeting  the  first 
two  criteria.  By  incorporating  in  them  certain  modifications  which 
will  be  described  later,  automatic  procedures  satisfying  criterion 
(ill)  are  derived  for  dividing  a disc  of  arbitrary  shape  into  quadri- 
lateral subdomains. 


I 

I 


I 


1 


i: 
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3.2.  The  Zienklewlcz-Phllllps  Technique 

A full  treatment  of  the  Zlenkiewicz-Phillips  mesh  generation 
method  Is  available  elsewhere  [5,10-12].  We  apply  the  technique  here 
to  decompose  an  arbitrary  simply-connected  disc  whose  boundary  is 
determined  by  a finite  number  of  points.  The  idea  is  to  divide  the 
disc  into  a number  of  zones  each  of  which  is  defined  by  eight  vertices. 
Each  of  these  zones  is  then  subdivided  into  quadrilateral  domains  in 
such  a way  that  the  zones  are  connected  through  those  quadrilateral 
vertices  (or  nodes)  which  lie  on  the  inter-zone  boundaries. 

3.2.1.  A Procedure  for  Subdividing  a Single  Zone 

The  square  [-1,  1]  x [-1,  1]  in  the  C-H  plane  is  mapped  invertibly 
onto  the  quadrilateral  domain  ABCD  in  the  x-y  plane  bounded  by  four 
parabolic  areas  on  which  eight  points  lie  (Figure  1).  These  eight 
points  are  the  images  of  the  points  in  the  ^-n  plane  with  coordinates 
(-1,  1),  (0,  -1),  (1,  -1),  (1,  0),  (1,  1),  (0,  1),  (-1,  1)  and  (-1,  0), 
in  that  order.  The  points  2,  4,  6,  8 in  the  x-y  plane  need  not  be  the 
midpoints  of  the  sides  on  which  they  lie,  but  must  lie  within  a certain 
neighborhood  of  the  midpoint  [12].  In  the  mapping,  functions  L^(5,ri), 

1 * 1,  2,...,  8 are  defined  such  that  the  position  vector  of  any  point  P, 
within  or  on  the  piecewise  parabolic  boundary  ABCD  can  be  expressed  in 
the  form 

8 

r (x,y)  - I L (C,n)  r (x,y)  , 

I'  i=l  ^ 


where  rj^(x,y)  - ir^. 


1.3 (x,y)  = r^,  etc. 
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The  function  L^(C.n)  are  such  that  n^)  = 6^^ , i=l,2,...,8; 

where  6^^  = 1,  i=j 

= 0,  1/j  (l.e.  the  Kronecker  delta). 

They  are  given  by  [10]: 


Lj^(C,n)  - 

- \ (i-C)(i-n)(i+C+n); 

L2(C.n)  = \ (i-n)(i-S^); 

LjCS.n)  - 

- i (IH)  (i-n)  (i-C+n) ; 

n)  = \ (i+C)(iV); 

LgCC,!!)  = 

- i (i+C) (i+n) (i-e-n) ; 

Lg(5,n)  =\  (i+n)(i-5^); 

L^C^.n)  “ 

- i (i-C) (i+n) (i+C-n) ; 

Lg(5,n)  = \ (1-5) (i-n^)  • 

(6a) 

By  writing  ^ = -1  + , s = 1,2, .. . ,2n,+l; 

S 1 

1 ^ t = l,2,...,2n„+l;  (6b) 

t n2  t- 

the  square  A’  B'  C D’  is  effectively  divided  into  4nj^  n2  rectangular 
subdomains  such  that  the  arcs  AB  and  CD  in  the  x-y  plane  are,  respec- 
tively, divided  into  2nj^  and  2n2  segments.  The  Z-P  mapping  guarantees, 
given  the  condition  mentioned  earlier,  that  the  image  of  a point  (s,t), 

with  coordinates  (5  ,0^).  in  the  C'O  plane  is  the  unique  point  P with 
s t st 

8 

8 

x(C^.  n^)  - n^)yi  • 


coordinates 


3.2.2.  Applications  of  the  Z-P  Procedure 


(1)  An  Arbitrary  Quadrilateral 

For  the  polygon  ABCD  (Figure  2),  taking  the  four  other  points 
required  In  the  Z-P  procedure  as  the  midpoints  of  the  sides  reduces 
the  parabolic  segment  of  Figure  1 to  straight  lines  segments,  and 
results,  on  applying  (6),  In  the  decomposition 

4 

Ip  = I M (C,n)  r,.  . 

1=1  ^ 


where  r^^^  = r^,  ^2  = r^,  = I^,  = r^,  with 


n)  = |-  (i-C)(i-n)  ; n)  = |-  (i+C)(i-n); 


n)  = |-  (1+0  (i+n);  m^(^,  n)  = I (i-0(i+n)  , 


and 


Mi(Cj,  Hj)  - i,j  = 1,2, 3, 4; 

r ^ = r (Cj^,  Tij^)  = iC-l.  1).  etc. 


(11)  A General  Polygon 

Consider  a general  polygon  with  N vertices.  If  N Is  even  then 
the  polygon  can  be  divided  into  p zones,  where  P “ J (N-2).  If  N Is 
odd,  then  by  introducing  an  extra  vertex  at  the  midpoint  of  the 
largest  side  of  the  polygon,  the  number  of  zones  is  y (N-1) . A 
polygon  with  p zones  then  has  n'vertlces',  where  n=N  for  N even,  and 
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n ” N+1  (N  odd),  and  is  represented  in  the  £;-n  plane  by  an  in-line 
arrangement  of  p squares  (Figure  3) . 

In  the  x-y  plane,  zone  r has  as  vertices  the  points?^, 

P^_^,  r ••  1,2,...,  p,  where  1'^  is  the  image  of  the  point  in  the 
C-0  plane.  By  dividing  P'^  ^r+1  ^^^to  2n^  segments  and  P^^j^  P|^  ^ 
into  2n2  segments,  the  entire  polygon  in  the  x-y  plane  is  reduced 
to  an  array  of  quadrilateral  subdomains  determined  by  (2nj^p+l)  x 
(2n2'hl)  vertices  as  nodes. 

3.3.  The  Gordon-Hall  Procedure 

The  Gordon-Hall  (or  G-H)  procedure  can  be  applied  to  domains 
whose  boundaries  consist  of  a finite  number  of  analytically-defined 
segments  through  the  use  of  several  blending  functions  [6,14,15]. 
Recently  the  blending  function  method  of  interpolation  has  been 
extended  to  the  semi  infinite  strip  [0,  1]  x [0,  “)  [19].  In  the 
present  work  we  restrict  our  application  of  the  method  to  finite 
simply-connected  domains  whose  boundaries  consist  of  at  most  four 
analytically-defined  arcs. 

In  the  G-H  procedure  the  unit  square  A'  B'  C'  D'  in  the  s-t 
plane  (l.e.  in  the  domain  S)  is  to  be  mapped  unlvalently  onto  the 
closed  domain  R in  the  x-y  plane  bounded  by  the  closed  curve  ABCD 
so  that  A,  B,  C,  D are,  respectively,  the  Images  of  A',  B',  C' , 
and  D*.  The  technique  assumes  the  existence  of  a continuous  vector- 
valued function  ^(s,t)  which  maps  the  s-t  plane  onto  the  x-y  plane 
such  that  the  boundary  A'  B'  C'  D*  is  mapped  onto  the  curve  passing 
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through  A,  B,  C,  D.  A function  U Is  then  constructed  which  maps  S 
onto  the  closed  region  R in  such  a way  that  this  function  Is  Identical 
to  F^(8,t)  on  the  boundary  of  the  square.  This  is  due  by  the  introduc- 
tion of  four  scalar  blending  functions 


<(>^(8),  ipj(t),  1 - 0,1;  j = 0,  1 , defined  by 


(j)^(s)  » 1 - s. 

<t)^(s)  = s , 

0 £ s £ 1 ; 

l^^(t)  - 1 - 1, 

= t , 

0 < t £ 1. 

In  taking  s =0,  s,  =1,  t =0,  t =1,  we  have  the  relations 
o 1 o o 

r l,k  = 0,  1; 


L j.  £ = 0.  1. 


Next,  projectors  P , P are  defined  as  follows: 
s c 

P^[F]  = <l>„(8)  F(8  ,t)  + (t)  (s)  F(s  , t)  , 

8 — O — O 1 


Pj!]  ' 1(8,  t^)  + F(s,  t^)  . 

The  required  function  lJ(s,t)  is  then  given,  in  terms  of  the  Boolean 

sum  P « P of  P and  P by  the  relation 
St  s t ■' 

U(s,t)  - (P„  • PJ  m - P„[F]  + P [F]  - P P.[F1, 

o w 8 t SC 

where  the  product  projection  P P is  defined  by 

8 t 


(7) 


P P III  ■ I I F(8.,  t.) 

® i-0  1-0  ^ J ^ J 


Whilst  the  right-hand  side  of  (8)  interpolates  to  £ at  the  four  vertices 
of  the  square  [0,  1]  x [0,  1]  in  the  s-t  plane,  the  vector  valued  function 
lJ(s,t)  given  by  (7)  interpolates  to  £ at  all  points  on  the  boundary  of 
the  square.  At  the  same  time  it  maps  the  points  with  coordinates  (s,t)  in 
the  interior  of  the  square  onto  the  nodes  with  position  vector  lJ(s,t)  in 
the  interior  of  R.  If  the  parametric  representations  F(s,t)  of  the  arcs 
AB,  BC,  CD  and  DA  in  the  x-y  plane  are  known,  then  (7)  yields  the  re- 
quired mapping  of  the  square  onto  the  region  R in  the  form 


U(s,t)  = (1-S)  F(0,t)  + s F(l,t)  + (1-t)  F(s,0)  + t F(s,l) 


- (l-s)(l-t)  F(0,0)  - (l-s)t  F(0,1)  - s(l-t)  F(1,0)  - st  F(l-l). 


3.4.  Mesh-Grading  Features 


3.4.1.  The  Gordon-Hall  Square 


In  Figure  5 the  solid  lines  divide  the  Gordon-Hall  square  into 
4nj^n2  rectangular  subdomains  with  Nj^(-  2nj^)  and  N2(“2n2)  subdomains 
along  A'B'  and  A'D',  respectively.  The  centers  of  these  subdomains 
lie  on  the  dotted  lines  as  shown. 

Figure  6 shows  a typical  subdomain  with  its  center  which  has 
coordinates  (s^^,  t^),  0 £ s^,  t^  £ 1,  1-1,2,...,  j-1,2,...,  N2. 


I 


The  subdomain  has  sides  2a^,  2a^ , Figure  6.  Using  the  notation 

+1  * ^1’  ^n  " ^n  +1  ~ ^2’  introduce  section  ratios  [20] 
and  where  a’  = a^  = k^A^,  in  such  a way  that 

the  following  linear  relations  hold. 

Ni  N 

a^^  = Aj^  - a(-2 — i)»  i “ ; 

N N 

= A^  - B(i  - - 1),  i = -^  + 1,,..,N^; 

where  a and  6 are  constants  to  be  determined.  Similar  relations 
hold  for  the  a j . On  expressing  a and  S in  terms  of  the  section 
ratio  kj^,  we  find  that 


2(1  - k )(^  - i) 

^i  " ^1^^ 


= A^[l  - 


2(1  - k^d  - -f 


1=1  2 — 

L 2 


-]  , i = + 1 N, 


r 1 

where  the  relation  I a,  • ^ immediately  gives 
1=1  ^ ^ 


“l  (1  + kj)Nj  • 

When  kj^  = 1 (i.e.  for  the  case  of  a uniform  subdivision  along  the  s-axis) , 
we  have  A^^  = 2^»  as  expected.  For  the  rth  Interval  along  the  x-axls, 
the  left-  and  right-end  points  are,  respectively,  at 
r-1 

\ 2a^  and  S2^=  + 2a^,  when  a^  is  given  by  (10).  The  center 

of  this  Interval  is  at  ^ ®r'  applying  a similar  procedure 


I. 
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to  intervals  along  the  t-axis  it  can  be  readily  verified  that  we 
obtain  a subdivision  of  the  square  in  the  s-t  plane  such  that: 

(i)  the  coordinates  (a.,  t.)  of  the  center  P'  of  the 


'i’ 


subdomains  are  given  by 


Si  = N— r-2  f2i(N^k^  - 2 + i - k^i)  - + 2],  i-1,2,..,  , 

1 ^1 

= 2 2(h^:2)  [(21-N^)(3N^  - N^k^^  - 4k^  - 2 i - 2ik^)  + 4kj^-2Nj^] 


i = ^ + . 


where 


Si  “ k^,  Aj^,  i)  , i = 1,2,...,N^,  say, 


tj  - SjCN^,  k2,  ’ j “ 1.2 N2. 


(1  + k2)N2  ’ 


(11a) 


(11b) 


(ii)  the  coordinates  of  the  vertices  of  the  subdomains  are 


(5j^.  Oj)  where 


2A  (i-1)  N 

^i  “ - 2 ^^1*^1  " 2 + i - kj^i),  i - 1,2,...,  -j- 

A (21  - N - 2) 

“ 2 ^ — 2-(^2)-  - N^k^  - 4k^  - 2i  + 2ik^), 


i - + 1,...,  + 1, 


^i  " ^1  ^**1’  '^l'  ^1’  ^ “ 1.2.-**.Nj  + 1 . 


(12a) 


(12b) 


(Nj.  k2,  A2,  J),  j - 1,2 N2  + 1 


(12c) 
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The  points  and  the  vertices  n^)  give  rise  to  images  in 

the  ..-y  plane,  under  the  mapping  U,  so  that  the  boundary  of  the  domain  R 
is  approximated  by  a polygon  R'  whose  vertices  are  the  images  of  the 


vertices  Hj)  which  lies  on  the  boundary  of  the  square.  This  polygon 

is  then  effectively  divided  into  quadrilaterals  by  connecting  the  image 


vertices  so  that  the  subdomain  containing  has  as  vertices  the  images 
of  the  vertices 


^ * ^^i+1*  ^ j ^ * ^^i+1*  snd  ^i+1^* 


i = 1,2,. j=l,2,...,N2. 


3.4.2.  The  Zienkiewicz-Philllps  Square 


Here,  since  the  Gordon-Hall  square  fO,  1]  {0,  1]  is  replaced 


by  the  square  [-1,  1]  x [-1,  1],  which  is  also  subdivided  into  rec- 
tangular domains  of  the  type  shown  in  Figure  6,  we  have  in  place 


of  (11)  and  (12); 


(i)  The  Coordinates  (s^,  t^)  of 

A N 

^i  “ “ ^ [2^(N^k^  - 2 + i - k^i)  - Nj^k^  + 2],  1*1,2 -f 

- 2(N^  " N^)(3N^-N^k^-4k^-2i  + 2ik^)  + 4k^  - 2N^], 


where 


1 (l+k^)N^ 


i = -^  + 1, . . . , 


(13a) 
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l.e. 


^1’  ’ 


(13b) 


tj  ■ ‘'2’  ^2’ 


(14) 


where 


A„  = 


2 (1  + k2)N2 


(ii)  The  Coordinates  (Cj^.  H^)  of  the  Vertices  of  the  Subdomains 

2A  (i-1)  N 

Cl  = - 1 + (Vl  - 2 + i - ^li),  1=1,2 


A (21  - N - 2) 

2(N.  --2)  »“l  - "l"!  - ‘‘‘l  - “ * . 


(15a) 


1 = -^  + I,-.-,  N^+1, 


or 


Cl  = Ci(Ni,  A^,  1)  ; 


(15b) 


rij  = 5j(N2,  k2,  A2,  j)  (16) 

By  setting  k^^  and  k2  to  values  less  than  1,  the  boundary  quadri- 
lateral subdomains  can  be  made  much  smaller  than  those  near  the  center 
of  the  domain  R' , so  that  the  charge  density  can  be  computed  at  points 
very  close  to  the  boundary,  without  the  need  for  the  very  large  values 


of  Ni  and  N2  which  would  be  required  with  uniform  subdivisions. 


I 

I 
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This  gives 


x(s,t)  = (1-s)  (l-t)x2^  + s(l-t)x2  + stx^  + t(l-s)x^ 

(17) 

y(s,t)  = (1-s)  (l-t)yj^  + s(l-t)y2  + sty^  + t(l-s)y^  . 

If  the  quadrilateral  is  a rectangle  with  sides  AB  = a,  BC  = b, 
and  with  its  vertices  at  A(0,0),  B(a,0),  C(a,b),  D(0,b),  then  (17) 
becomes 

x(s,t)  = sa,  y(s,t)  = tb  , 


so  that  the  Gordon-Hall  mapping  reduces  to  a two-way  stretch  of  the 
square  [0,  1]  x [0,  1]  parallel  to  the  x-  and  y-axes. 


3.5.2.  'Ihe  triangle 


Figure  8 shows  a triangle  with  vertices  A(xj^,yj^),  B(x2,y2)« 
C(x2.y2),  such  that  BC  is  the  longest  side.  Taking  ABC  as  a 
degenerate  quadrilateral  A^^  B^^  with  vertices  A^(x|,  » 

Bi(x^,  yp,  C^(xy  yp,  y^).  where 

x^  - Xj^,  X2  “ X2,  *3  = J (^2  ” *3’ 


yi.  y2 


■T^ 


Fig.  8.  An  image  triangle  in  the  x-y  plane 


Fig.  9.  A subdivision  of  the  triangle  in 
Figure  8. 


x(s,t)  = (1-s)  8(2-t)x2  ^ \ 


y(s,t)  = (1-s)  (l-t)yj^  + ^ s(2-t)y2  + y t(2-s)y^  . 

An  obvious  simplification  results  if  A is  taken  as  the  origin.  In  (18) 
we  observe  that  the  parametric  lines  s - constant  and  t * constant  are 
mapped  onto  straight  lines  in  the  triangle  ABC  (Figure  9) . 


3.5.3.  The  Ellipse 


The  points  A(a,0),  B(0,b),  C(-a,0)  and  D(0,-b),  which  are  the 

2 2 2 2 2 2 

extremities  of  the  axes  of  the  ellipse  bx  +ay  =ab  (Figure  10), 
are  chosen  as  the  Images  of  the  vertices  of  the  Gordon-Hall  square. 
Parametrically,  the  equation  of  the  ellipse  is,  for  0 < s,  t < 1: 


Along  AB:  F(s,t)  = F(s,r), 


SIT  , . STT 

X = a cos  — , y = b sin  — ; 


Along  BC: 


Along  CD: 


Along  DA: 


X » a cos  Y (t+l)i  y = b sin  j (t+1); 


TT  TT 

X = a cos  (s+1),  y “ - b sin  (s+1); 


X = a cos 


, y = - b sin 


Using  (9)  we  find 


t I 


! r 


U(s,t)  - (1-s) 


a cos  — (t+1) 


+ (1-t) 


-b  sin  ^ b sin  ^ (t“l) 


u ^ STT 

b sin  — 


giving 


a cos  j (s+1) 


-b  sin  2 (s+1) 


- (1-s) (1-t) 


(l-s)t  - s(l-t) 


7 («’t) 
a 


(1-s)  cos  — 


s sin  — + (1-t)  cos  — 


t sin  + 3+t-l 


y , TTt  /I  \ j ■'ft  . V . TTS  ITS  , 

t-  (s,t)  = s cos  (1-s)  sin  — + (1-t)  sin  t cos  ---  + t-s 

b I I 12 

Here  the  parametric  lines  in  the  s-t  plane  are  mapped,  in  general,  onto 
curved  lines  in  the  x-y  plane,  unlike  for  the  quadrilateral. 


3.5.4.  The  Ellipse  Sector 


If  the  angle  of  the  sector  of  an  ellipse  is  2a,  we  introduce  a 
fourth  vertex  at  the  front  C(a  cos  a,  b sin  a),  where  2a,  2b  are  the 
axes  of  the  ellipse.  By  defining  the  parameters  s and  t for  the  seg- 
ments AB,  BC,  CD  and  DA  of  the  boundary  as  shown  (Figure  11),  and 
observing  that  the  vector  valued  function  F(s,t)  assumes  the  following 
forms: 


AB:  F(8,0) 


F(s,l) 


we  obtain 


~a  cos  (ta) 

; 

F(l,t)  = 

_b  sin  (ta)_ 

a cos(2-s)a'l  at  cos  2a 

, DA  = F(0,t)  = 

b sin(2-s)a  bt  sin  2a 


— (s,t)  = s cos  (ta)  + t cos(2-s)a  - st  cos  a > 
a 


J (s,t)  = s sin  (ta)  + t sin(2-s)a  - st  cos  a . 

D 


3.5.5.  The  Circle  with  Two  Tangents  (’The  Circle-With-Tangents') 

Figure  12  illustrates  a disc  bounded  partly  by  a circle  of  radius  a 
and  by  the  two  tangents  drawn  to  it  from  a point  A,  the  tangents  being 
Inclined  at  an  angle  a to  the  x-axls.  These  tangents  meet  the  circle  at 
B and  D,  and  we  take  as  the  'corner'  nodes  the  points  A(a  cosec  a,  0), 

B(a  sin  a,  a cos  a),  C(-a,0),  D(a  sin  a,  - a cos  a).  The  parametric 
forms  of  the  various  segments  of  the  boundary  are  easily  seen  to  be 
(cosec  a + s(sln  a - cosec  a) 

AB:  £(s,0)  = a 

s cos  a 


BC:  F(l,t)  = a 


CD:  F(s,l)  - a 


DA:  F(0,t)  ••  a 


cos  { (j  - a)  + t (^  + a) ) 
sin  { (^  - a)  + t (^  + a)} 


cos  { (j  - a)  + s (^  + a) } 
-sin  { (-j  - a)  + s (j  + a) } 


cosec  a + t(sln  o - cosec  a) 
- t cos  a 


From  these  descriptions,  the  image  £(s,t)  of  the  point  (s,t) 


is  given  by 


(s,t)  = (1-st)  cosec  a + (s+t-2st) (sin  a.  - cosec  a) 


+ s cos  { (j  - a)  + t (y  + a) } + t cos  f (^  - a)  + s (j  + a) } 


- (s+t-2st)  sin  a + St; 


(s,t)  = s sin  { (y  - a)  + t(Y  + a) } 


7T  If 

t sin{a(-2  - a)  + s(j  + a)} 


+ (s-t)  cos  a. 


3.5.6.  The  Annulus  Sector 


A sector  of  an  annulus  is  considered  here  as  an  example  of  a 
non-convex  planar  domain.  Such  a sector  is  shovm  in  Figure  13,  the 
radii  being  and  R2  (R^^  < R2)"  The  parametric  equations  of  the 
boundary  can  be  written  down  immediately  in  the  form 


F(s,t) 


IRj^d-s)  + s R2}  cos  at 


{Rj^(I-s)  + s R^}  cos  at  J , 


where  a is  the  sector  angle. 

On  applying  (7)  to  (22)  we  find  that 


P (F)  » F(s,t),  whilst  P (F)  - P P (F)  = 0. 
s — — t~  st~~ 
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Thus  U(s,t)  is  identical  to  lF(s,t).  The  required  image  in  the  x-y  plane  is, 
therefore,  given  by 


x(s,t)  = {Rj^(I-s)  + s R2}  cos  at 
y(8,t)  ■■  {Rj^(I-s)  + s R2}  sin  at  . 


(23) 


Gordon  and  Hall's  result  for  the  quarter  annulus  [6]  is  easily  seen 
to  be  a particular  case  of  (23) . 

Figures  14(i),  to  14(vli)  show  typical  subdivision  arrangements 
for  the  case  in  which  the  G-H  square  is  divided  uniformly. 

3.5.7.  Contours  Defined  by  a Set  of  Points 

The  application  of  the  Z-P  method  to  a contour  defined  by  four 
points  was  given  in  detail  in  Sections  3.2.1  and  3. 2. 2.1.  For  contours 
with  more  than  four  defining  points,  which  must,  therefore,  be  divided 
initially  into  zones  (Section  3. 2. 2. 2),  the  final  decomposition  is  ob- 
tained for  the  first  zone  by  applying  the  procedure  for  a single  zone. 

For  zone  r(r  > 1)  the  process  is  the  same  except  that  in  using  the 
four  vertices  P',  P'.,  P'  ,,  P',,  defining  this  zone  (Figure  3)  the 
parameters  (C  , n,.)  given  by  (6b)  are  required  only  for  s-2,3, . . . ,2n,+l, 

SC  X 

t*l,2, . . . ,2n2+l,  since  the  inter-zone  nodes  corresponding  to  s-l  have 

already  been  accounted  for  as  those  corresponding  to  s ■ 2nj^+l  in  zone  r-1. 

Since  all  the  squares  are  in  line,  general  'coordinates'  (i,j)  can  be 


9 
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Fig.  14 (v).  A subdivision  scheme  for  a sector  of  an 
annulus . 


Fig.  14 (vi).  A sector  of  an  annulus:  a subdivision 

scheme  with  some  matchpolnts  lying  out- 
side their  subdomains. 


m. 


T- 

\ 


i 


35 


adopted  to  identify  the  vertices  for  the  p zones,  where  i=l,2, , . .2nj^p+l, 

J*l,2, . . . ,2n2+l. 

The  layout  of  the  nodes  and  subdomains  generated  in  a regular 
hexagon,  for  = 6 = and  = 1.0  = k^,  are  shown  in  Figure  lA(vill). 


4.  THE  INTERACTION  MATRIX  AND  THE  POTENTIAL  INTEGRAL 


4.1.  Integration  over  Triangular  and  Quadrilateral  Domains 


The  rectangle  shown  in  Figure  6 will  be  mapped  onto  a curvilinear 

quadrilateral  subdomain  of  the  disc  in  the  x-y  plane  (Figure  15)  so  that 

the  centre  of  the  rectangle  is  mapped  onto  a point  P^j  within  the 

quadrilateral.  This  quadrilateral  has  as  vertices  the  points  Q^j t 
2 3 4 

^ij  ’ ^ij  ’ ^ij  ’ are,  respectively  the  Images  of  the,  points  in  the 

original  square  with  coordinates  Dj),  ’^j+1^ 

and  (5^,  Hj).  For  integration  purposes  the  four  points  in  the  x-y  plane 

are  taken  as  the  vertices  of  a polygon,  so  that  the  disc  is  ultimately 

replaced  by  a collection  of  such  quadrilaterals  (Figure  14) . 

Using  the  notation  of  Section  3,  Figure  15  may  be  taken  as 

illustrating  the  rth  subdomain  (2^  of  the  disc.  The  four  triangular 

subdomains  k=l,2 4,  are  those  obtained  by  joining  P^^  to  each 

of  the  vertices  of  the  quadrilateral,  so  that  has  as  vertices  the 
k k+1  5 1 

points  Q^j,  , P^j,  where  . It  has  been  shown  that  the 

piecewise  integration  over  indicated  in  (5a)  depends  on  Integral 
expressions  of  the  form  These  expressions  are  best  evaluated  by 

the  use  of  the  area  (or  barycentric  ) coordinates  of  the  triangle 


Fig.  15.  A quadrilateral  subdomain  in  the  x-y  plane  divided 
into  four  triangles. 


A 


Fig.  16.  Area  coordinates  for  the  triangle  ABC. 
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The  use  of  these  coordinates  is  a standard  procedure  in  the  literature  [22-25], 
These  coordinates  are  usually  denoted  byC,  H,  C.  OjiCt  H,  C£l.  and  are 
defined,  for  any  point  P(C»n,5)  in  the  triangle  ABC  (Figure  16),  by 


where  A.__  Is  the  area  of  triangle  ABC,  etc.  Thus  ^ + n + C = 1,  and  the 
ABC 

equations  of  the  sides  AB,  BC  and  CA  can  be  expressed,  respectively,  in 
the  forms:  C = 0,  C“0,  ri  = 0. 

For  any  function  f(x,y)  which  is  integrable  over  triangle  ABC, 

1 1-n 


f(x,y)dxdy  “ 2 A. 


f {x(C,n),  y(C,n)}  dCdn, 


(24a) 


where  denotes  the  domain  of  the  triangle. 

The  transformation  , il'  = H converts  (24a)  to 


f (x,y)dxdy 


d-n’)  f{C(C',  n’).  n'}  dC'dn' 


(24b) 


This  transformation  effectively  maps  the  triangle  bounded  by  ^ - 0,  n “ 0, 
5 +ri+l  » 0 in  the  ^-n  plane  onto  the  square  [0,  1]  x [0,  1]  in  the  C'-n' 


plane. 


It  is  easily  shown  that  if  the  vertices  of  the  triangle  are 


A(Xj^,  yj^),  B(x2,  yj).  C(x^,  y^)  then 
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— 

- 

X 

Xi  X2  X3 

y 

“ 

yi  y2  y3 

n 

1 

111 

C 

^ — 

— — 

l.e.  X = 

y “ + b^n  + 

where  - x^,  ” ^3*  *^3  “ *3» 

The  inverse  relation  Is 


(25) 


”5' 

^2  " ^3  "'b  ■ ^^2  ’'2^3  “ V3~ 

1 

n 

S — 

2A 

^3  ~ ^1  ’^l  - *3  ’^3^1  " ^^1^3 

y 

c 

..  . 

yi  y2  X2  - X3  x^y2  - y^X2 

1 

where  2A  - (x^^  - X3)(y2  “ 3^3)  “ (*2  “ ’'3^^^!  " ^3^*  I^l 
the  area  of  the  triangle. 

Equation  (26)  can  be  expressed  in  the  form 


(26) 


5 - Aj^x  + Bj^y  + ,13“  A^x  + B^y  + C2  , 


C - 1 - ? - n, 

where  A^^,  B^^,  etc.  are  constants. 

For  quadrilaterals  of  the  type  shown  in  Figure  15,  uniqueness  of 
the  area  coordinate  representation  of  points  in  the  various  triangles  can 
be  ensured  across  the  four  interior  boundaries  by  arranging  that  such 
lines  should  correspond  to  either  C “ 0 or  n “ 0 in  each  of  the  two 
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triangles  in  which  It  belongs.  If  the  coordinates  of  a typical  point 


In  the  triangle  are  rij^,  then  the  lines  can  be 


represented  as  * 0 * n^,  whilst  has  the  representation 


“ 0 ■ etc.,  as  in  Figure  17. 


4.2.  The  Potential  Due  to  a Uniform  Charge  Distribution 


over  a Triangular  Domain 


The  potential  at  a point  (x',  y')  in  the  plane  of  the  triangle 

s s 


f2  . , due  to  a unit  charge  density  on  this  triangle  is  (4ire  )~  C , , 
rx  o srK 

where  is  given  by  (5b).  On  applying  the  transformation  leading 

to  (24b),  and  dropping  the  primes  on  ^ and  n»  we  have  (appendix. 


Section  A2.1.1.) 


1 1 


■ “rt  I I y(g.n))-  “rk  I ‘K-  <”> 


0 0 


where  A is  the  area  of  £2  , and  g(x',  y';n)  is  given  by 


“rk 


■ log  {a+bn  +V(a+bn)^  + (c+dn)^}  - log  {a*+b'n  + V (a'+b’n)^(c+dn)^} 

(28) 


where  a^j^-  aj^^  + a^^^  ; a - 1 + R(x;,  y');  b - P-1; 


c - S(s;,  y;);  d - Q;  a'  - R(x;,  y;);  b'  - P; 


^ “^®lrk’’lrk  ®2rk’’2rk^^“rk’  ^ ” ^®lrk'’2rk  " ®2rk*’lrk^^“rk 


R(x;,  y;)  - - x;)  + a2^^(c2^^  - y;)}/^  ^ 


S(x;,  y;)  - {aj^^(c2^^  - y^)  - ’ 


i r 


i I- 


Fig.  17.  A quadrilateral  subdomain  in  the  x-y  plane,  with  area 
coordinate  representation  of  interior  and  exterior 
boundary  lines. 


Brk 


Fig.  18.  A triangular  subdomain 
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Here,  on  referring  to  Figure  18: 


^Irk  ~ ’'irk  ” *3rk’  '’irk  " ’'2rk  " ’'srk*  ‘’irk  ’‘Srk’ 


^2rk  “ ^Irk  " ^3rk’  '’2rk  " ^2rk  " ^3rk’  ‘’2rk  “ ^3rk 


The  evaluation  of  the  n-integral  (Section  A2.1.2)  gives 


a . 

(2!^)  Cgrk  " * 

rk 

where  E(ri)  = E(a,b;n),  E'(n)  = E(a',  b';n)  , 

and  E(a,b;ri)  = ^ log  {a+bn  + V (a+bn)^  + (c+dn)^} 


ad  - be 

~2  2 
<k/b  + d 


log 


d{a+brH- V(a+bn)^-r(c+dn)^}-  (c-Mn)  (b4  Vb^-M^) 


d{  a+brH" 


V(a+bn)' 


+(c+dn)''}  - (c+dnXb-Vb^+d^) 


An  alternative  expression  for  which  is  less  susceptible  to 

round-off  errors  is  (Section  A2.1.2) 
a , 

(2!^)  C^rk  ° *'■  ~ ’ 

rk 


(29) 


where  F(n)  = F(a,b;n),  F'(n)  = F(a',b';n)  , 


and  F(a,b;r|)  = log  {a+bn  + V (a+b  )^  + (c+dn)^) 

a 


+ log  [ (b^+d^)n  + ab+cd  + V { (b^+d^)n+ab+cd}^+(ad-bc)^] 

dV  b^+d^ 


For  an  arbitrary  field  point  with  coordinates  (*g»yg»^a)»  which 
need  not  lie  in  the  plane  of  the  triangle,  the  distance  function  D in 


(27)  is  replaced  by 
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5 {x’  y’,  z’;  x(C,n).  y(C.n)}  = {(x-x')^  + (y-y')^  + 


8 8 


S 8 


In  place  of  (28)  we  now  have  (Section  A2.2.1) 


“tk  K- 


log  (a+bri  + V (a+br|)^+(c+dr|)^+e^}  - log{a'+b'n+  V(a'+b'ri)^+c+dn)^+e^} , 


where  e = z'/ot  , and  the  other  constants  are  as  for  (28).  On  evaluating 
s irR 

the  n-lntegral  (Section  A2.2.2)  we  find  that 


rk 


s(K>  yl»  z!;n)cih  = g(i)  + g'(0)  - g(0)  - g'(i), 

s s s 


where  G(n)  = G(a,b;n),  G'(n)  = G(a',b';n), 


and  G(a,b;ri)  * ^Qg  {a+bn  + V(a+bri)^  + (c+dn)^  + e^} 

a 


+ ~ log  Kb^+d^in  + ab+cd+V((b^+d^)n+ab+cd)^+(ad-bc)^+e^<b^-M^)l 

dV^ 


. e . -1 
+ -T  sin 
a 


de)/(a+bn)^  + (c4dri)^  + e^ 


[y{(c+dn)^+e^}  {ad-bc)^+e^(b^+d^)}J 


Thus  the  potential  at  an  arbitrary  point  In  space,  due  to  a unit  charge 
density  on  the  triangle  Is  where 


I 


u 


1 [ 
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rk 


This  expression  for  reduces,  as  expected,  to  (29)  when  e “ 0 

J 


(30) 


+Vx^  + Y"]  occur  In  both  (29) 


Expressions  of  the  form  log  [X 
and  (30).  If  X is  negative  such  expressions  are  best  computed  by  using 
the  alternative  form 

(— 

,.2 


log 


- X + y x^ + 


RESULTS  AND  OBSERVATIONS 


5.1.  Charge  Distribution 


The  charge  density  on  an  ellipse  with  semiaxes  a and  b is  [21] 

O'  = (1  - ^ - ^)  . (31) 


_5_  (1 . 

4llab  2 .2' 

a b 


where  (x,y)  are  the  coordinates  of  a typical  point  in  the  ellipse, 
and  Q Is  the  total  charge.  If  an  elliptical  disc  is  maintained  at 
unit  potential,  then  observing  that  the  charge  accumulates  on  both 
sides  of  the  disc,  (31)  gives  the  charge  distribution  a on  the  disc  as 


. _ z!)-1/2 

2tTab  2 ,2^ 

a b 


(1 


(32) 


r 

it 


where  C Is  the  capacitance  of  the  disc  as  given  in  the  appendix 
(Section  A.I.). 


I The  exact  charge  distribution  (32)  has  been  compared  with  the 

computed  values  for  a circle  of  radlus/cm  (Figure  19),  and  for  an 
i ellipse  with  a “ 5 cm,  b = 1 cm  (Figure  20).  The  results  show  that 

j for  a given  number  of  subdivisions  of  the  ellipse  or  circle  (l.e.  for 

fixed  values  of  and  N2  as  defined  for  the  Gordon-Hall  square)  the 
j use  of  non-uniform  section  ratios  and  gives  a charge  distribution 

which  is  more  accurate  than  for  uniform  section  ratios.  As  expected,  a 
! non-uniform  section  ratio  allows  charge  densities  to  be  computed  very 

near  the  edges,  and  the  approach  to  the  square  root  singularity  at  the 
edges  is  more  clearly  brought  out.  For  = 8 = for  example,  at 
comparable  points  in  a circular  disc,  working  with  uniform  section 
ratios  (l.e.  k^^  = 1.0  = k2)  results  In  error  margins  of  24%  near  the 
1 edges  and  1.7%  near  the  center,  whereas  for  k^  = 0.1  = k2  the  error 

j margins  are  4%  and  0.7%.  For  an  elliptical  disc  with  a = 5 cm, 

b = 1 cm,  the  figures  are  52%  and  1.4%  for  unity  section  ratios,  and 
j 4%  and  0.8%  for  section  ratios  of  0.1. 

The  distribution  for  unity  section  ratios  for  a square  of  edge 

i 2 cm  (Figure  21)  agrees  with  Harrington's  result  [3]  to  within  1% 

j everywhere  except  at  the  node  nearest  to  the  edge  where  the  difference 

is  1.7%  (Figure  21).  The  beneficial  effect  of  non-uniform  section 
j ratios  is  again  evident.  In  Figure  22,  drawn  for  a regular  hexagon  of 

side  1 cm,  the  distribution  shows  the  expected  symmetry  and  the  rise 
} in  the  charge  density  as  one  approaches  the  edges. 

I Figure  23  compares  the  computed  charge  distribution  over  a 

semicircle  with  that  for  a full  circle.  A similar  figure  (Fig.  24) 

^ I 

4, 

1 J * 


CHARGE  DENSITY  (X  lO"* COULOMBS / M 
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has  been  drawn  for  a square  and  a triangular  half  of  this  square. 


5.2.  Capacitance 

It  is  shown  in  the  appendix  that  the  capacitance  of  an  elliptical 
disc  with  semi-major  axis  a is 


C = 4iTae  /K(e)  , 
o 


(33) 


where  K is  the  complete  elliptic  Integral  of  the  first  kind  and  the 
modulus  e is  the  eccentricity  of  the  ellipse.  A computer  program 
(PROGAE)  computes  C from  this  expression,  using  the  Hastings  approxi- 
mation [17,28]  for  K(k)  in  the  form 


where 


K(k)  = I a k^’'  - 2[  I b k^*"]  log  k + e(k), 
r=0  ^ r=0  ^ 

|e(k)|  <2  X 10~®,  k^  = 1 - k^  , and 


Sq  » 1.38629436112 

= 0.09666344259 

a^  “ 0.03590092383 

- 0.03742563713 

a,  - 0.01451196212 
4 


b^  = 0.12498593597 

b^  = 0.06880248576 

b^  = 0.03328355346 

b,  = 0.00441787012 
4 


For  a circle  of  radius  a we  use  the  well-known  result  C = 8ae  , 

o 

which,  as  shown  in  Section  A.l,  is  a particular  case  of  (33). 

Convergence  curves  which  show  the  variation  of  C with  (N^^N^)  ^ 


for  the  circle  and  the  square  (Figure  25)  indicate  that  section  ratios 


SQUARE  (2  CM  X 2 CM) 
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of  0.1  give  a more  accurate  value  for  the  capacitance  than  unity 
section  ratios.  The  greater  accuracy  obtained  with  section  ratios 
less  than  unity,  which  was  observed  for  the  charge  density,  is  also 
reflected  in  all  the  capacitance  computations.  This  Improvement  in 
accuracy  agrees  with  similar  observations  made  in  the  investigation 
of  current  distribution  effects  in  perfectly-conducting  cylinders  [20]. 

To  estimate  the  capacitance  for  an  infinite  number  of  sub- 
divisions (l.e.  the  asymptotic  capacitance)  we  use  a power  series  pro- 
cedure which  gives  C as 

" C 

C - I ^ (34) 

1=0  (N^N^) 

If  there  are  m sets  of  computed  data  in  the  form  {N,  , N„  ; C }, 

Is  zs  s 

8 » then  by  truncating  (34)  in  the  form 


m-1 

C=  I 

L= 

the  following  matrix  equation 


1=0  (N^N^) 


1 ’ 


r-l 


'Vl'  ■ ■ 


r = 1,2,...,  m;  s = 1,2,...,  m, 

is  then  solved  to  give  as  an  estimate  of  the  asymptotic  capacitance. 

For  the  circle,  the  asymptotic  capacitance  computed  for  various 
sets  of  data  corresponding  to  Figure  25  are  given  in  Table  1.  The 


deviations  of  the  values  for  section  ratios  of  1.0  and  0.1  are, 
respectively,  vdthln  0.9%  and  0.2%  of  the  exact  value.  The  capaci- 
tance for  = 10  = N2  (l.e.  for  the  largest  values  of  and 
used)  Is  found  to  differ  from  the  exact  value  by  0.37%,  whereas  for 
“ 8 « N^and  for  the  same  section  ratios,  the  error  is  slightly 
less  (0.31%).  This  gives  an  indication  of  the  effect  of  the  increase 
of  round-off  with  the  number  of  subareas. 

Since  an  inscribed  polygon  underestimates  the  area  of  any 
convex  domain,  one  would  expect  the  computed  capacitance  for  such 
domains  to  be  lower  than  the  exact  value.  This  Is  found  to  be  the 
case  for  the  circle  and  for  the  other  convex  domains  investigated. 

Tables  2,  9 and  12  list  the  computed  values  of  capacitance  for  the 
ellipse,  the  semicircle  and  the  'circle-with-tangents'  and  also 
give  the  percentage  deviation  of  the  area  of  the  inscribed  polygon 
(computed  by  adding  the  areas  of  all  the  triangular  subdomains)  from  the 
area  of  the  original  domain. 

For  polygonal  discs  on  the  other  hand,  the  mesh  generation 
procedures  ensure  that  the  inscribed  polygon  is  the  same  as  the 
original  domain.  Figs.  14(1),  14(vll)  and  14(viil),  so  that  the  area 
of  the  inscribed  polygon  as  computed  agrees  with  the  area  of  the 
original  domain  up  to  the  fifth  significant  figure  at  least  (see 
Tables  4-7  for  the  regular  hexagon,  the  diamond,  the  rectangle  and 
the  triangle).  Computations  for  the  square  of  edge  2 cm  (Figure  25 
and  Table  3)  show  that  the  capacitance  for  = 10  ” N2  and  section 

ratios  of  0.1  (l.e.  0.8102  vy  F)  is  within  0.43%  of  the  asymptotic 


I 

I 
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value  (0.8137  UW  F).  The  asymptotic  value  obtained  by  Table  4 is 
0.34  e.s.u.  which  is  0.8156  up  F.  For  the  diamond  (Table  5a)  and 
the  rectangle  with  dimensions  3.3  cm  ^ 1.4  cm  (Table  6),  the  capacitance 


computed  for  = 8 = N2  and  0.1  section  ratios  are,  respectively, 

0.6510  PU  F and  0.9029  pp  F.  De  Meulenaere  and  Van  Bladel  [29]  give 
0.661  pp  F and  0.897  pp  F. 

The  results  for  the  sector  of  an  annulus  (the  only  non-convex 
domain  Investigated  here)  indicate  that  care  must  be  taken  in  approxi- 
mating domains  by  polygons.  It  is  important  to  ensure  that  the  match 
point  generated  lies  not  only  within  the  curvilinear  quadrilateral 
corresponding  to  it,  but  also  within  the  final  quadrilateral  obtained 
when  the  curved  sides  are  replaced  by  straight  lines.  When  the  sector 
angle  is  180“  and  the  section  ratios  are  unity,  this  condition  holds 
when  “ 8 “ N^,  but  it  is  violated  when  = 4 = N2  (Figure  14  (vi)). 
A low  section  ratio  can  also  give  a match  point  lying  outside  its 
quadrilateral.  Thus,  for  example,  for  0 = 180“  and  = 8 = N2 
(Table  10)  a section  ratio  of  0.1  produces  this  anomaly,  as  also 
happens  when  0 = 120“  and  = 4 = N^.  Match  points  are,  however, 
properly  situated  in  the  latter  case  when  the  section  ratios  are  unity. 

Generally  a violation  of  the  match  point  condition  at  certain 
points  leads  to  a distortion  of  the  charge  distribution,  including 
the  existence  of  negative  charge  distributions  at  some  of  these  points. 
It  also  gives  for  the  approximating  polygon  an  area  which  is  larger 
than  that  of  the  original  domain  mainly  because  of  the  resulting 


1 

I 

I 


duplication  of  triangular  subdomains  as  the  area  is  computed.  As  a 
guide,  well-positioned  match  points  can  be  ensured  by  choosing 
and  in  proportion  to  the  lengths  of  the  corresponding  arcs  of  the 

domain  boundary. 

The  effect  of  the  distortion  in  the  shapes  of  the  quadrilaterals, 
when  the  match  point  condition  is  met,  has  been  investigated  by  comput- 
ing the  percentage  errors  in  the  values  obtained  for  the  area  and  the 
capacitance  for  ellipses  with  aspect  ratios  (i.e.  minor  axis/major  axis) 
from  1 to  0.2  (Table  11).  The  results  show,  for  the  coarse  subdivision 
used  (N^  = 4 = N^),  that  the  area  is  in  error  by  2.6%  in  all  cases.  The 
corresponding  figure  for  the  capacitance  varies  from  5.4%  for  an  aspect 
ratio  of  1 to  4.9%  when  this  ratio  is  0.2.  These  errors  are  thus  due 
more  to  the  coarseness  of  the  subdivision  than  to  the  shapes  of  the 
quadrilateral  subdomains. 

The  computations  also  show  that  the  capacitance  of  half  of  a disc 
is  greater  than  half  the  capacitance  of  the  complete  disc. 

6.  CONCLUSION 

A method  of  moments  technique  of  solving  the  singular  Integral 
equation  satisfied  by  the  charge  distribution  on  perfectly-conducting 
flat  plates  of  arbitrary  shape,  which  uses  automatic  mesh  generation 
techniques,  has  been  presented.  The  method  relies  on  the  use  of  graded 
quadrilateral  subdomains  and  has  been  applied  to  a variety  of  disc  shapes. 
Including  the  rectangle  for  which  the  quadrilateral  subdomains  reduce  to 
rectangles.  The  procedure  yields  charge  distributions  and  values  of 


I 

1 


57 


i 


f 

capacitance  which  are  In  excellent  aRreemont  with  the  exact  values  i 

for  the  circle  and  the  ellipse,  as  well  as  with  published  numerical  ■ 

results  for  the  rectangle  obtained  by  ocher  methods  of  computation. 

The  results  Indicate  that  with  a moderate  number  of  subdomains  It  is 
possible  to  obtain  values  of  capacitance  which  are  accurate  to  within 
0.5%. 


, ' 

i 


Table  1; 

The  Circle 

Radius  = 1 cm;  Capacitance 


Exact  Value  of  Capacitance:  0.70800  uijF 


- N2  <•  N 

■ 

■=  k2  = k 

Asymptotic  Capacitance 

Computed  Value  (PVF) 

% Error 

6,  8,  10 

1.0 

0.7017 

0.89 

6,  8,  10 

0.1 

0.7035 

0.65 

6,  8 

0.1 

0.7065 

0.21 

Table  2; 

The  Ellipse 

Semi  Axes:  3 cm.  1 cm:  Capacitance 


0.15307 

0.15450 

0.15607 

0.15495 


% Error 


Capacitance 

Computed  Value  (upF) 

% Error 

1.753 

4.9 

1.817 

1.5 

1.802 

2.3 

1.838 

0.29 

Asymptotic  Capacitance 


Computed  Value  (UPF) 


% Error 


1.818 

1.846 


Table  7 
The  Triangle 

Ad.lacent  Sides:  2 cm,  1 cm;  Included  Angle;  90** 
Exact  Area:  0.20000  x ]0 


N 

k 

Area  (x  10 

Capacitance  (yuF) 

A 

1.0 

0.20000 

0.56A8 

A 

0.1 

0.20000 

0.5958 

6 

1.0 

0.20000 

0.5800 

6 

0.1 

0.20000 

0.6050 

8 

1.0 

0.20000 

0.5879 

8 

0.1 

0.20000 

0.6077 

Asymptotic  Capacitance 


N 

k 

Computed  Value  (iJUF) 

A,  6.  8 

1.0 

0.6001 

A,  6,  8 

0.1 

0.6109 

Table  8 

The  Semicircle:  Radius  1 cm 

-3  2 

Exact  Area:  0.15708  x 10  m 


N 

k 

Area  (x  10 

Capacitance  (MhF) 

A 

1.0 

0.15307 

0.A855 

A 

0.1 

0.15616 

0.5105 

8 

1.0 

0.15607 

0.50A1 

8 

0.1 

0.15537 

0.5186 

Table  9 
The  Quadrant 
Radius  = 1 cm 

-4 

Exact  Area:  0.78540  x 10  m 


Area  (x  10  m ) j Capacitance  (HMF) 


0.78413 

0.78273 


0.3507 

0.3612 


Table  10 

The  Sector  of  an  Annulus:  Radii:  R^  = 1 cm,  R 


(a)  Sector  Angle 
Exact  Area: 


1 Area  I 

-3  2 

Computed  (x  10  m ) 

% Error 

0.45922 

2.6 

0.45897 

2.6 

0.53833 

14.2 

Capacitance  (yuF) 


0.9456 

0.9502 

0.9399 


(b)  N = 


Area 


% Error 

Capacitance  (uyF) 

5.5 

0.90000 

4.5 

0.7113 

2.6 

0.5994 

Table  11 

The  Ellipse:  Semi-Axes 


N - 4.  k = 1.0.  b 


Variation  of  Capacitance  with  a/b 




Area  (x  10 

tn^) 

Capacitance  (UPF) 

Computed 

Exact 

% Error 

Computed 

Exact 

% Error 

0.30615 

0.31416 

2.6 

0.6695 

0.70800 

5.4 

0.36737 

0.37699 

2.6 

0.7350 

0.77719 

5.4 

0.42.860 

0.43982 

2.6 

0.7980 

0.84365 

5.4 

0.45922 

0.47124 

2.6 

0.8288 

0.87604 

5.4 

0.48983 

0.50265 

2.6 

0.8591 

0.90794 

5.4 

0.55106 

0.56549 

2.6 

0.9185 

0.9703 

5.3 

0.61229 

0.62832 

2.6 

0.9766 

1.0314 

5.4 

1.5307 

1.5708 

2.6 

1.753 

1.8436 

4.9 

Table 

12 

The  Circle-with-TanRents: 

Tangent  Angle  (a' 

Radius  (cm) 


Computed 

Exact 

0.37320 

0.38202 

0.38264 

0.38264 

1.5210 

1.5144 

1.5306 

1.5306 

Capacitance  QJy  F) 


0.7509 

0.7869 
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APPENDIX 


A-1.  The  Capacitance  of  an  Ellipsoid 

The  capacitance  C of  an  ellipsoid  with  semlaxes  a,  b,  c Is  given 

. o ,.21 

by  Smythe  as 


Sire 


dx 


(x  + a^)(x  + b^)(x  + c^) 


(i 


In  order  to  evaluate  the  Integral  In  (1)  we  use  the  following  result 

26 

given  by  Jahnke,  Emde  and  Losch  : 

00 

- ^ 0 

/(b^  + t^)(c^  + t^) 

1 Hi 7 

where  cot  (J)  ■ x/c,  k*  — yc  -b, 
and  F(<ti,k)  “ 


du 


i ! 2 T 

0 yi  - k sin 


the  Incomplete  elliptic  Integral  of  the  first  kind.  Under  the  trans- 
2 2 

formation  t - x + c , the  integral  In  (Al)  becomes 


I 


af ^ 

J t,  1 ■ ' "2 2772 72 27 

cY(t  + a - c)(t  +b  - c) 


Using  (2),  and  assuming  that  a > b > c,  we  obtain  Innediately 


68 


The  capacitance  of  the  ellipsoid  can,  therefore,  be  expressed  in  the 
form 


- .^,2  2.1/2.  , ^ -1.,  2 -2,1/2,  2 ,_2,l/2,  2 2-1/2 

C = 4'rTe^(a  -c  ) /F[sln  (1  - c a ) . (a  - b ) ' (a  - c ) 


This  expression  corrects  the  result  given  by  Smythe 


21 


(A3) 


The  Ellipse 

When  c = 0,  because  F(Tf/2,k)  = K(k)  (i.e.,  the  complete  elliptic 
integral  of  the  first  kind  with  modulus  k) , (A3)  becomes 

C = ATrae^/K(e),  (A4) 

where  a is  the  semi-major  axis  of  the  ellipse  and  e the  eccent 

For  a circle  of  radius  a,  putting  a = b in  (A4)  gives  imF-f  d-:  it . ly 

C = 8ae  , 
o 

a well-known  result. 


A, 2 The  Potential  due  to  a Uniform  Charge  Distribution 
Over  a Triangle. 

We  present  here  a summary  of  the  procedures  adopted  in  deriving 
the  relations  (27)  and  (30)  for  the  potential  due  to  a uniform  charge 
distribution  over  a triangle.  The  notation  adhered  to  is  that  used  in 
Section  4.2,  unless  otherwise  indicated. 


I 


A. 2.1  Field  points  In  the  plane  of  the  triangle. 


A. 2. 1.1  The  g- Integral. 


We  have 


g(x'  y';n)  = 

S o 


(1  - n)dg 


5 D{x^,y^;x(g,n),y(C,n)} 


where 


D{x^,y^;x(g,Ti)p^(C,n)}  = Y{x(g,n)  - xM  + {y(g,n)  •-  y^}  , 
x{g,n)  - ” ^1^  ®1’  ~ ^2^  ■*■  ®2' 


” ®lrk^^  ~ ®1  “ ’’irk'^  '^Irk’  ‘^Irk  ‘ ‘^Irk  " ’^s 


A2  » a2rk(l  - n).  Bj  = b2^1^n  + ^2rk  “ ‘'2rk  " 


We  then  have 


D^{x' ,y’ ,x(g,Ti)  ,y(5.n)}  “ (A^  + A^){(g  + p')‘  + q‘} 
s s i.  ^ 


2 . .2, 


.^2  . _2, 


where 


“rkd  ” ^>P'  ' ^^^Irkhrk  "■  ^2rk'’2rk^  ^rk^lrk  ®2rk‘=2rk’ 
°rkd  " ’^>‘1  “ - a2rk^lrk^  ®lrk‘'2rk  " ®2rk‘^lrk* 


These  relations  give 


“rk8d;.y;;n) 


0 t/<C  + p’)^  + 


“rk^di.YgJn)  “ [cosh 


■ii±je:  / 


(A-5) 


I r'. 


I 

I 


which  Is  the  same  as  (28) 
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A.2.1.2.  The  n- Integrals 


In  order  to  obtain  (29)  we  require  a primitive  for 
f(n)  ■ log{a  + brj  + y (a+brj)^  + (c+dji)^}  - log{a'+b'n  + ]/(a'+b'ri)^  + (c  + dn)^} 

It  Is  easy  to  see  that  such  a primitive  has  an  algebraic  form  which  does 
not  depend  on  the  sign  of  (c  -H  dp).  Thus,  provided  - c/d  ^ [0,1],  we  may 
consider.  Instead,  a primitive  for  the  function 

f'(n)  - + 1)  - + 1) 

However,  It  can  be  shown  that  even  when  - c/d  e [0,1],  the  Improper 

Integral  I f'(il)  dn  converges.  It  Is,  therefore,  sufficient  to  obtain 
^0 

a primitive  E^(a,b:ri)  lor 

for  the  case  c + dn  > 0. 

On  writing  a + bn  “ (c  + dn)sh  6,  we  find  that 


I f(a,b; 


_v  ad  - be  r 6 I _d0__i 

“ d ^dshO  - b " I dshe  - b^ 


and  It  Is  easily  shown  that 


d6 

dsh0  - b 


2 2 
b + d 


dCshe  + cosh  9)  - (b  + V b^  + d ) 
d(sh6  + cosh  6)  - (b  - V b^  + d^) 


so  that 


Ej^(a,b,n)  - E(a,b;n)  +-^-5-^  log  (c  + dn). 


A primitive  for  f(a',b';n)  Is 


E^(a'  ,b'  ;n)  ■ E(a'  ,b'  ;n)  + log  (c  + dn) 
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Thus  the  required  primitive  for  f(n)  is 

E(n)  “ E(a,b;n)  - E(a',b';n).  (A6) 

This  primitive  yields  the  first  expression  for  given  in  Section  4.2. 

Alternative  Expressions 

We  now  consider  alternative  methods  of  simplifying  the  expression 


H(a,b;n)  “ log 


d{a  + br)  Vca  + bn)  -HCc+dri)  } - (c -Hdr))  (b  4-Vb^ +d^) 


I I 2 7~  r5 — T" 

d{(a  + bn)+  VCa  + bri)  +(c  + dr|)  } - (c  + dri)  (b -Vb  +d  ) 

(i)  Direct  reduction 

Now 


H(a,b;n)  “ log 


d(cosh  e + sh  6)  - (b  + Vb^  + d^) 

2 


d(cosh  6 + sh  6)  - (b  - Vb^  + d 


Let 


••  d cosh  0 +'\/b  + d , 


••  d sh  6 + b. 


Then 


Thus,  H(a,b;n)  “ log 


• d cosh  6 - Y b + d , 


Yi  + B^ 

Xi+^1 

- log 

i.e.  H(a,b;n)  “ - log 


2 2 
1*1  - 6? 


6, 


“ d sh  e - b. 


(Yi  - a^)  (Y^  + Bj^) 
(Y^  - otj^XX^  + Bj^) 


+ constant. 


On  substituting  for  X^  and  B^  and  ignoring  constants  in  the  primitive 
from  now  on,  we  find  that 

T— 7 , 


H(a,b;n)  ” - log 


12  2 

d -<•  b ah  6 + cosh  6 Vb  d 
d sh  0 - b 


(A7) 


1 


9^ 


Since  (c  + dn)(dsh0  - b)  ■ ad  - be,  a constant,  we  find,  on  simplifying 
the  right-hand  side  of  (A7)  that 


T 
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2 2 

H(a,b;ii)  = - log[(b  +d  )ri  + ab  + cd  + 


V 


[ (b^  + d^ri  + ab  +cd}^  + (ad  - bc)^  ] 

(A8) 


(li)  The  use  of  subsidiary  Integrals. 


We  obtain  here  two  primitives  hj^(n)  and  h2(r|)  for  the  function 


I 2 1 

h(a,b;n)  = -y  (a  + bn)  + (c  + dn)  , 


(A9} 


where 


and 


hj(n)  - j (c  + + 1 dn 


h . -i/h^  + f (n  + ab  + cd.^  ra  + c (ab  + cd)  ■, 

h2(n)  ih  + d J in  + 2 2^  ^ 2 . ,2  “ ..2  . ,2.2^ 

■'  b=d  b+d  (b+d) 

For  hj^(n),  putting  a + bn  = (c  + dn)sh  6,  as  before,  gives 


hj^(n) 


s — f (b  + d )n  + ab  + cd}  t/ (a  + bn)^  + (c  + dn)^ 

2(b  + d ) 


(ad  - be) 

I(b^+d^r^^ 


log 


d{a  + bn  + V 

^(a +bn)^  + (c +dn)^}  “ (c  + dn)  (b  + V 

d{a  + bn  + V 

^a  + bn)^+ (c  + dn)^}  - (c+dn)(b  -t/ 

When  the  integral  for  h2(n)  is  calculated,  however,  we  have 


h„(n)  ^ j { (b^  + d^)n  + ab  + cd)  '^(a  + bn)^  + (c  + dn)^ 

2(b  + d ) 


■'  logl(b^+d^)n  + ab  + cd+  ^{(b^  + d^n  + ab+cd}^  + (ad-bc)^] 

2(b  +d  )^' 

Theae  two  expressions  for*  the  .Rriniitlye  of  h(a,b;n)  thus  establish  the 
equivalence  already  implied  by  (A7)  and  (A8). 


A. 2. 2 Arbitrary  field  points 


A. 2. 2.1  The  t- Integral 


Here,  Instead  of  (A5) , we  have 

1 


(1  - n)  dS 


Dfx’  y'  z';x(£:,n),y(?,n)} 

s s s 


where,  on  proceeding  as  In  Section  A. 2. 1.1  we  find 


rK  s s s 


'(C  + p’)  + q 


-2  2 . e 

where  q “ q + 

(1  - n) 


.22  ,2. 

and  a . e = z' 

2 rk  s 


.2  . 2, 


“rk®^*s’ys’*s’'^^  “ log{a  + bn  + V (a  + bn)  + (c+dn)  + e } 


-log{a'  + b'n  + Y(a'+b'n)^+(c  + dn)^  + e'^}. 


.2  . 2i 


with  a,  b,  a',  b',  c,  d,  being  defined  as  In  Section  A. 2. 1.1 


A. 2. 2. 2 The  n- Integral 

(1)  The  p-n  transformation  and  applications 
We  define  a new  variable  p as  follows: 

p “ -/(c  + dn)^  + . 

This  allows  us  to  proceed  as  In  Section  A. 2. 1.2,  so  that  we  need  only 
consider  the  primitive  of  the  function 


I • (a .bin)  - log  { - + 1 } . 


The  relation 


a + bn  ■ p sh  ((> 
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gives 


r2  7 

pdsh4!  = k + bYp  - e, 


k = ad  - be, 


so  that  p satisfies  the  x-quadratlc 


2 2 2 2 7 ? ? 

X (d  sh  (j)  - b ) - 2xdksh  (J)  + be  + k =0. 


We  choose  the  root  p given  by 


/~T  2 

up  = dk  sh  4)  + b Y k - e u, 


2 2 2 . 

where  u = d sh  4i  - b , since  this  expression  for  p reduces  to  c H 
for  field  points  which  lie  in  the  plane  of  the  triangle  (i.e.  for  < 
The  primitive  for  f'(a,b;ri)  is  thus 


G(a,b;n)  = 4)n  - 

Since  bn  = -a  + p sh  ()),  we  have 


n d(j). 


G(a,b;n)  = G - £ g + G,  , 
Obi  I 


where 


and 


Let 


2 

G = 1-P  sh  <|)  = i + ^ f 

^’0  b • ^’1  d ^ d J 

G2  = I d<t. 


d4i 

~1  2 T 

d sh  d)  - b 


r~2  I I — 2 — Z 

X=dsh<|>+V^>  + Y = dsh((!-yb  + d, 

a”dsh(|)  + b,  B = dsh4)  - h 


Then  XY  = aB,  and  we  find  that  on  writing 

2 


-G 


• 1 1 ' 


1 d 


d4 


'1  I .2  ,.2  ^2 

d sh  4(  - b 


2bVb^  + 


Y + B , Y+  a , 
[log  Trr^  - log  — — ], 


x + 6 


X + a 


(A-10) 

dn 

= 0). 


(All) 


ai 


uyiCk  -e  u)(b  +d  +u)  k/b  +d  |kVb+d  +u- 


These  results  ultimately  show  that 


G(a,b;n)  = 


+ dri)<l 

d 


. 2 ,2 
b + d 


- log 


I 
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(11)  Alternative  expressions 

The  expression  L = - L2  + occurs  In  (A13),  where 

T 1 Y + g 

° X + a ’ 

, _ , Y + 3 

h ° X + 3 ’ 


12  2 (2  2 

kd  cosh  <p  + Vb  + d yk  - eu 


kd  cosh  (p  - vb 


2 ^ ,2  /72  2“ 

+ d vk  - e u 


(A-16) 


We  give  here  a simplified  form  of  L which  is  similar  to  the  expression 
H(a,b;ri)  defined  by  (A8). 


(a)  Direct  reduction 


Let  y = V (a  + b ) +(c  + d)  + e. 


It  is  easy  to  write  (A16)  in  the  alternative  form 


where 


L3  = + L22. 


p + d^  (b  cosh  (p  + yb^  + d^  shi}))  - kd 

r~2  2 /~2  ? 

p Vb  + d (b  cosh  (p  - yb  + d sh()))  + kd 


(A17) 


L32  - log 


,2 

Vb  + d 


d sh  - b cosh  0 


12  2 

b + d sh  ij)  + b cosh  (}> 


-A-furtbM  reduction  gives 

L31  *•  p(a,b;Ti)  - log  {k^  + e^(b^  + d^)},  

2 2 n 2 

where  P(a,b;ri)  " 2 log  [ (b  + d )r)  + ab  + cd  + y yb  + d ], 

Li  . .2  lr..2  . .2.  . . . .i2  . , . ..^2  . .2,^2  . ,2, 


(A18) 


. ■ t 


and 


y 


Using  the  relations  (All)  it  is  easy  to  show  that 


- . Xg - Yg 

^32  Xa  - Y6  ’ 


Li  + L2  = log  - ; 


. . , , . , , Vb  + d - 2b 

and  that  L,  - = log  , a constant. 

1 2 /~2  I 

Vb  + d + 2b 


(A19) 


(A20) 


We  thus  have  immediately 


- L2  + = P(a,b;r))  + constant, 

so  that  the  function  + L^,  being  a component  of  the  primitive 

G(a,b;ri)  given  by  (A13),  is  equivalent  to  P(a,h;n).  The  function 
P(a,b;ri)  is  very  similar  to  H(a,b;n),  and  reduces  to  it  when  e = 0. 


Reduced  expressions  for  and 

The  relations  (A17)  - (A21)  imply  that,  to  within  arbitrary 
constants , 

2 2 2 2 
1 , a - Y . 1 , B - Y 

L = -T  log  — =■  , L = - log  —y , and 

^ ^ X^  - a ^ ^ X - B'^ 


(A21) 


13=  P 


(a,b;n)  + log  /b^  + sh  l b_cosh 


)/b^  -b  d^ 


d sh  + b cosh 


On  substituting  for  X,  Y,  a,  B,  <(>  and  P,  the  following  explicit  reduced 


forms  result : 


“ log 


b(a  + bn)  - d /(c  + dn)^  ^ + y Vb^  + d^ 


(A22a) 


d(a  + bn) 


+ b J 


(c  + dn)  + e 


.11 


L2  = log 


d(a  + bri)  - b V(c  + dn)^  + ^ 


b(a  + br|)  + dV(c  + dn)^  + + y ■\/i 


, 1. — Z 
y Vb  + d 
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(A22b) 


and 


= 2 log[(b^  + d^)  + ab  + cd  + y l/b^  + d^] 


+ log 


2 2 2 2 2 
b { (c  + dr|)  + e } - d (a  + bri) 


(2b  +d^)(a+bri)  + b {(c  + dr))^  + e } + 2b (a  + br))yl/b^  + d^ 

(A22c) 

(b)  Auxiliary  Integrals 

In  line  with  the  procedure  based  on  the  function  h(a,b;ri)  defined 
by  (A9),  we  consider  here  two  primitives  for  the  function 


J(a,b;ri)  = -^  (a  + bn)^  + (c  + dn)^  + e^. 


where  the  primitives  Jj^(n)  and  J2(n)  are  given  by 


Ji(n)  = 
J,(n)  “'s/i 


P ^ ■'■1  dn ; 


2 2 
b +d 


I 2 2 2 2 2 

/r_  . ab  + cd-i  , ra  + c (ab  + cd)  , e ■, 


With  the  expression  for  p given  by  (AlO),  we  find  that 


■^1^^^  ” ^2b  ^(a  + bn)^  + (c  + cn)^  + 


/l2  .2  2,  f 

^ 2b J 


- f 

2d  J 


de  , .,2 

A A "T  Dk 

d sh  0 - b 
du 


J 

2 

e_ 

2d 


de 

2 ^ 2 2 
(d  sh"  e - b ) 


du 


/(b^ +d^  + u)(k^  - e^u)  ""  ^ u l/(b^  + d^  + u)  (k^  - e^u) 


: hi 


I 


It  Is  easy  to  show  that 


(dsh  0 + b)  a(b  + d ) (b  + d ) 


(dsh  0 - b) 


d cosh  0 
2 

B(b  + d^ 


(6  + d )■ 


These  relations,  together  with  (A12),  give  the  primitive 

j (n)  = Uc  + dn)(b^  + d^)  +__kbj.  7(3  + bn)^  + (c  + dn)^  + 
^ 2b(b  + d ) 

. {k^  + e^b^  + d^)}  ,,  TXT! 


However,  on  evaluating  the  integral  for  J2(n)f  flrid 


. ((c-^dnXb^^dV.kb)}  V(a  * bn>'  - <C  + dn)'  -b 


2b  (b^  + d^) 


{k^  4-  e^(b^  + d^)}  , . 


Thus  L - L„  + is  equivalent  to  P(a,b;n),  as  before. 


i 


A. 2. 2. 3 Planar  Interactions 


The  primitive  G(a,b;r|)  given  by  (Al3)  can  be  written  in  the  form 


1.  ._\  _ (c  + dn)4>  . e , -1  , ed  cosh  (t 

G(a,b;Ti)  = t + -y  sin  [ 1 

d d — — 

Vk  + e (b  + d ) 


[L^  - 1^2  + L^]. 


2d  Vb  + d 


For  planar  interactions,  e = 0,  so  that  X = Y = a=  6 = 
and  <|)  = 0.  Also, 

Xi  Xj  + 6^ 

L,  - log  ^ , and  L,  + L-  = log  — r , so  that 
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Lj^  - L2  + Lj  = ~2L2 , a relation  which  can  also  be  established  using  the 
explicit  forms  given  by  (A22). 

Thus,  G(a,b;n)  = E^(a,b;n),  showing  that  when  the  field  point  lies 
in  the  plane  of  the  uniformly  charged  triangle,  we  recover  the  primitive 
for  planar  interactions  given  by  (A6). 


